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Abstract — We consider the problem of operator identification 
in quantum control. The free Hamiltonian and the dipole 
moment are searched such that a given target state is reached 
at a given time. A local existence result is obtained. As a by- 
product, our works reveals necessary conditions on the laser 
field to make the identification feasible. In the last part of this 
work, some Newton algorithms are proposed together with a 
continuation method to compute effectively these operators. 

I. Introduction 

In the last decades, quantum control has known significant 
improvements both at theoretical and practical levels (cf.[l], 
[2], [3], [4] and references therein). Results have been 
obtained on existence of controls [5], [6], [7], [8], [9] or effi- 
cient ways to compute and carry out laser fields that achieve 
some goals concerning the state of quantum systems [10], 
[11], [12], [13], [14], [15]. On the other hand, the design of 
relevant laser fields plays also a major role when the goal 
is to identify some properties of the quantum system to be 
controlled. In this way, some methods have been designed 
to identify finite dimensional systems characteristics [16], or 
to compute discriminant laser fields [17]. 

Note that operator identification in relation to the 
Schrodinger equation has already been studied in the litera- 
ture. As an example, we refer to [18] for a theoretical result, 
where no laser interaction is considered. 

In this paper, we focus on the case where only one 
fixed laser is used to identify in finite given time the free 
Hamiltonian and the dipole moment. From the theoretical 
point of view, we obtain a local existence result: we prove 
that the inversion is always possible in the neighborhood of 
some particular states. As a by-product, we emphasize some 
features of the laser fields that enables the identification. 

Following the local approach we use to obtain this result, 
we present in a second part, a time discretized setting and 
fixed-point methods to solve numerically our problem. In 
particular, a Newton method is proposed together with a 
continuation method that allows us to solve problems where 
the local assumption does not hold. 

This paper is organized as follows: the mathematical 
formulation of our problem is given in Section |II] and a 
local controllability result is presented in Section |III] In 
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Section IIVI we present the algorithms to solve numerically 
the identification problem. We conclude with some tests in 
Section |Vl 

Let us finally introduce some notations concerning par- 
ticular matrix sets that will be used throughout the paper 
Given Na £ N, we denote by C^'''^'' and R^-'^^-' the sets of 
matrices of size Nd x Nd with complex and real coefficients 
respectively. Then, define 



u = 


= {MeC^'^'^^ M*M^MM* 


s = 


= {Mec^<^'^^ M* = M}, 


Sr = 


= {MeM^'''^^ M* = Af}, 


cO _ 


= 5r n {Af e R"'""^ Mk.k^o, 



= l,...,Nd}, 

where M* denotes the adjoint matrix associated to AI and 
Id is the identity matrix of C^'*'^''. Here, for the sake of 
simplicity, we have omitted the dependence of these sets 
with respect to Nd- In what follows, we denote by Kz and 
Sz denote respectively the real and imaginary parts of a 
complex number z. Given a matrix M, we denote by Af^ 
its transposed. 

II. Setting of the problem 

Fix r > 0, and consider a system U{t) G U whose 
dynamics over [0, T] is ruled by the Schrodinger equation: 

iU{t) - [Ho + e{t)^]U{t), (1) 

(7(0) = U^mt, (2) 

where Hq e iSr is the matrix of the internal Hamiltonian, 
e(t) e i2(0,T;R) a laser field, fi £ Sm the matrix 
associated with the dipole moment. For relevant applications, 
the matrices Hq and /i are not supposed to commute. The 
initial state Uinit is fixed. In this equation, e is given and the 
pair {Ho, /i) £ iSr x 5^ is searched such that at time t — T, 
the state reaches a given target state Utarget, i-e., 

U{T) = Utarget. 0) 

In other words, given the mapping 

ip: 5r X 5g -^ U 

iHo,n) ^ UiT), 

the main question that will be investigated in this paper is 
the surjectivity of ip. 

In our work, the internal Hamiltonian H^ is searched as 
real Hermitian (i.e. symmetric) matrix. This is a particular 
situation as in general it is only supposed to be complex 
Hermitian and not real. Nevertheless, for the applications 
we have in mind this restriction is very natural since the 
Hamiltonian is a sum of a kinetic operator and a potential. 



both real. For the same reasons, we suppose that the dipole 
moment /i is real (Hermitian thus symmetric) but we assume 
moreover that the diagonal elements are null. This additional 
assumption is motivated both by invariance properties (the 
diagonal of Hq as matrix commutes with the diagonal of 
/i as matrix) but also by the desire to identify an unique 
pair {Ho,fi) since in this way the number of unknowns 
(dimension of 5r plus that of S^) equals the number of 
equations (the dimension of U). 

Note that one can easily prove the following conservation 
property: 

Vie[0,T], \\U{t)\\u = \\U,mt\\u, 



where we have denoted by 
scalar product 



the norm associated to the 



{A,B) eUxU^tr{A*B). 

This problem is related to inverse problems in quantum 
control [17], but unlike previous works, we do not aim here 
at designing relevant laser fields to identify the pair {Hq, fi) 
but rather to investigate the properties of the fields e{t) 
that make Equation Q invertible and algorithms to compute 
numerically the corresponding solution operators Ho and /i. 

III. Local controllability result 

In this section, we present some theoretical results about 
the local inversion of Equation (|3]l. More precisely, we make 
use of the calculus of variations to obtain a local inversion 
theorem. 

Given a pair {Hq, fi), we first introduce the tangent space 
•AHa,ij,, which is the space of matrices defined by: 

^ffo,M = {^^ e C^'''^^ M*U{T) + U{T)*M = 0} . 
We then consider the differential operator of ip defined by: 

d(p{Ho, n) : SrX Sm^ -^Ho.i^ 
{6Ho,6fi)>^dU{T), 

where SU{T) is solution at time t = T of the linearized 
Schrodinger equation: 

iSU{t) = [Ho + e{t)^i\5U{t) + [SH^ + e{t)5^i]U{t), 

and U{t) follows equation ([T]l. 

We will prove that Lp is an onto mapping using the fact 
that dip also satisfies this property. This strategy is motivated 
by the following known result: 

Theorem 1: Supposed that d(p{HQ, /i) is an onto mapping, 
i.e. 

W e AHo,t.,^iSHo,5fi),dipiHo,fi){dHQ,dfi) = V. 

Then ip is locally onto in a neighborhood of {Hq, /i). 
We shall prove that dp is an onto mapping on the neighbor- 
hood of all states of the form Uq :— p{Hq,0) G U. To do 
this, we compute explicitly an inverse mapping. 

Theorem 2: Given H^ G 5e, define Vb as the matrix that 
diagonalizes Uq := p{Ho,0) in the following way: 

Uo{t) = Fo*e^^(*-¥Vo, 



with A the diagonal matrix with coefficients Aa G M, a € 
Nd, 1 < a < Nd- Suppose that for a 7^ 6, 1 < a < iV^, 1 < 
b<Nd, 



Aa 7^ Afc 



(4) 




(i)e*«A„,.(t-i:)^^j _^Q^ (5^ 



Then dp{Ho, 0) is an onto mapping and its inverse is given 
by: 

ip : V e^Ho:M '^ {SHo,6fi). 

The matrices SH^ and Sjj, are given by: 

5Ho := Vo*SHoVq, S^i := V*6t^Vo, 

where the coefficients ha.b and nia^b of the matrices SHq 
and (5/i are given by: 



ma,b 



"^JVa^b 



-a,b 



^Va,b - 7r^'=iVa,b 



sin((5Aa,6 — 



J. S\a,b i/a / b ^^^ 



rria 
hn 



T 



if a = b. 



Here Va.b: ci,b ^ Nd, I < a < Nd are the coefficients of 

iVQ*UQ(T)*V'Vo and e^;,, := ^ (j^ e{t)e^^^-->''^*-i'>dty 
Proof: We fix V' G Aho.^ and solve 

dp{Ho,0){SHo,Sfi)^V'. (7) 



First, one can show the identities: 

piHo,fiydipiHo,0)iSHo,Sfi) = UoiT)*SUo{T) 

^~if Unity i6Ho+eit)S^i)Unit)dt, (8) 



where the variation SU^ is defined by the evolution equation: 

i6Uo{t) = [Ho + e{t)fi]SUo{t) + [SHo + e{t)Sti]Uo{t). (9) 

Note that such an identity holds also when /i ^ 0. Since 
Uo{T)* is invertible, showing that O has a solution is 
equivalent to show that 



Ua{ty{SHo+e{t)6fi)Ua{t)dt = V, 



(10) 



has a solution, with V := iUo{TyV' G S since V G 
AHo,ti- A nice property of the trajectory t H> Uo{t) is that 
Equation (fTOl i can be solved explicitly. Indeed, let us denote 
by Ua,6, ha^b and nia.b, with a, & G N, 1 < a, 6 < iV, the 
coefficients of the matrices VqVVq, VoSHqVq and Vq^/xVq* 



respectively. Expanding ( flOl l gives rise, in the case a ^ b to where i„ = ^42^(-ffo + SnfJ-) 



VaM - K., / e'(^-^^)(*-*)di 



+ma,6 / £(t)e'(^"-^'')(*~*)di 



T 
sin((5Aa,b — ) 

= ha^b FT Vma,b£{5\a,b), 

OAa,b 

where 6Xa,b = K - h and e{SXa,b) 

In the case a = b, one finds that 

T _ T f'^ 

Va,a = ha,b-^ + ma,b£(0) = ft-a.hlT + rila.b / e(t)dt. 



Let us now detail the effect of variations 5Hq, dfi in iJp and 
yit on the sequence {Un)n=Q Nt- We have: 

{Id + Ln)5Un+l + SLnUn+l = {Id - L„)6Un 

-SLnUn, 

SL^{U„+i + U^) = {Id-L^)SU^ 

-{Id + L,,)5Un+i, 

{Un+l+Unr5Ln{Un+l+Un) = -2{U*^^^5Un+l 



'U*JU,,), 



where 5Ln — ^^^^{511^ + EnSf-i). This finally gives rise to: 



Un+lSUn+l - U*5Un 



-I AT {dHo + £„(5/i) 



c/„ 



■ArpV^ {Un+l + Un)* ,^ ^ ,Un+l 

=-zAr2^ 2 {SHo + SnSfi) 

n=0 



Note that t he assum ption SHo, S^ £ Sm combined with Since the initial value is fixed, we obtain: 
£(^Aa,6) = e{SXb.a) impHes that Va,b — Vb.a, so that V £ S. 
The result follows. ■ 

Remark 1: In this theorem, we have defined rua.a arbi- 
trarily. 

This theorem gives a first hint about conditions required to 
identify (i/o, /^)- Condition (|4| is weaker to the standard non- 
degeneracy condition 

V(a,fe)^(a',6'), \b-Xa^K'-\a', 

and is in practice often satisfied. Condition (|5]l deals with 
the laser field itself. It is a non-resonant condition to control 
the system. 

IV. Numerical methods 

In this section, we present two algorithms to solve Q. The 
strategy we follow is a direct adaptation of previous results 
and proofs: we consider local approximations based on fixed 
point iterative solvers. In our approach, a crucial step consists 
in obtaining an appropriate time discretized version of dl). 
In the first part, we build such an approximation that enables 
the exact computation of the derivative of the final state 
U{T) with respect to (iJo,A*) and derive from this setting 
a numerical strategy. 

A. Time discretization 

In order to simulate numerically Equation ([T), we intro- 
duce the following time discretization: give Nt G N, we 
T 

denote by AT — the time step and for n = 0, • • • , Nt 

Nt 
by Un and e„ the approximations of U{nAT) and e{nAT). 

In order to preserve the unitary property of the matrices U{t) 
at the discrete level, we use a Crank-Nicholson scheme ruled 



Un 



(11) 



This result can be seen as a discretized version of (|8) where fi 
is not necessarily null. We insist on the fact that such a result 
is specific to the Crank-Nicholson discretization. As far as we 
know, no other numerical solvers give rise to discretization 
of dSj where the variations SHq and Sfi are explicit. 

B. Fixed points methods 

We now present some iterative solvers to compute solu- 
tions of Q. 

1) A Newton Method: In the discrete setting, we still 
denote by (p the operator: 



f 



{IIa,fi) f-> Unt 



To solve the equation ip{IIo, ^) = U tar get, a Newton method 
would consist in the following iteration: 



d^{Hl /•) • {SH^Jfi'^) = - {^{H^, /•) - Utarget) , 



(12) 



where k is the iteration index, SHq = Hq — H^, S^ 



,k+l 



In our case, ( fT2] i reads: 



^^ Nt ~ U target ~ I^ATt 



Using (fTTT l. one can rewrite this equation as follows: 



AT 



by the formula: 

. Un+l — Un 



Nt-1 ijTk , TjkX 






n=0 



{Ho + Sn^J.) 



Un+l + Un 



l{{U%^)*Utarget-Id), 



AT ' "^^ 2 

The corresponding iteration is then given by: 

{Id + Ln)Un+l = {Id - Ln)Un 



where we recall that the unknowns are SH^ and 5^^ . This 
equation has generally no solutions, since its left hand side 
belongs to S what is not the case for its right hand side. To 
solve this problem, we replace i {{U^ )* Utarget — Id) by 



a first order approximation 5''"' G iSr. Two possible choices 



exp(-zS''') := 



(U 



'U, 



target 



.{U^YU, 



target 



TT* TT^ 
'-^ target^ Ni 



(13) 
,(14) 



In the numerical tests, the same behavior is observed when 
using the first or the second definition. 

Remark 2: The previous method can be simplified to 
obtain a procedure where no matrix needs to be assembled 
and the inverted during iterations. Instead of up-dating at 
each iteration in the pair {Ho,fi) in the term dip{Ho,fi) 
of Formula ( fT2] l. one can keep a constant approximation 
{H^''^, ^f^'f) of the solution. We denote by (C/;'=/)„^o ... _^^ 
the corresponding sequence of states. The iteration then 
reads: 



N^ 



-1 (Tjref 

AT Y, +' 



u^^fy 



{5H^ + e^jfJ/i* 



[/, 



re/ 



Jjref 



^S' 



where S'' is defined in the previous section, see (fT3T l and (fT4l i. 
Note that such an algorithm is actually a time-discretized 
version of the fixed point used in the proof of Theorem |2] 
except that here // is not supposed to be null. 

2) Implementation of the iterative solvers: Both previous 
methods require inversions of linear systems which are not 
given explicitly in our formulations. To fill in this gap, we 
explain here how to assemble the matrices, i.e. to rewrite the 
equation 



AT 



JVt-1 (J J 
n=0 



(dHa + SnO^) 



-5, 



in terms of linear system. In what follows, we denote by 
Xm the vector representation of a matrix M consisting in 
concatenating vertically its columns. A first step to do this 
is to note that the later equation reads as follows: 



(Nt-1 
n=0 



/Nt-1 \ 

\XsH„+ATi E£„Mc/„^,,Jx,^ 



Xs, (15) 



with 

M, 



X fcron(C/J+i/2,lAr, 



Here, kron denotes the Kronecker product, Un+1/2 = 
"^2 " ' '^1'^ '■^'^^ ^y ^^^^ product of two matrices A and 
B is denoted hy A. x B and Ijv^ denotes the matrix of 
-^Nd.Nd \yjjose coefficients are equal to 1. 

A second step must then be carried out: since the matrices 
5Hq and 6fi are symmetric, one has to consider the columns 
of the matrices in (flST l that correspond to the coefficients 
of 6H0 located, e.g., above the diagonal and the coefficients 
of (5/i located strictly above the diagonal. In the same way, 
only the lines of the resulting system that correspond to 
the coefficients located above the diagonal of S shall be 
considered. 



Taking the real and the imaginary part of the equations, 
the resulting system is of size Nj. 

C. A continuation method for global controllability 

The algorithms proposed in Section IIV-BI are only locally 
convergent. The purpose of this section is to present a 
continuation method that enables to extend their range of 
application. 

As mentioned above, numerous methods exist to solve the 
control problem where the laser term e in Equation ([T]l is 
unknown and iJo and /i are given [15], [11], [10]. Based 
on this fact, the method we propose is the following. Given 
an initial guess (iJQ,/i*'), find a control e° such that U^^, 
the final state associated to [H^^ji^) correctly approximates 
Utarget- Givcn 6 G [0,1] , wc define the interpolated fields 
£* = (! — 9)e'^ + 9e. A fixed point method as the one 
presented in Section lTV-Bl can then be applied with {Hq, fjP) 
as an initial guess to solve the operator control problem with 
£^. Our algorithm consists in repeating this procedure by 
solving iteratively the operator control problem associated to 



the field e"^" using (Hi; 



-\/-i) 



as initial guess. Carrying 



this procedure up to = 1 enables to solve the original 
problem. 

V. Numerical results 

In this last section, we present some numerical results 
obtained with the algorithms of the previous sections. As 
a laser term in Equation ([T), we use e{t) = sin(t). The other 
numerical data are Nd = 5, To = 10, Nt = 10^, T = 27rro 
and AT = T/N. 

A. Newton Method 

We first test our Newton method. In this way, we choose 
randomly a pair {HQ,fi), with coefficients in [—1,1] and 
compute the corresponding final state Unt- Then, we start 
the Newton procedure with an initialization {Ho + AHo,iJ,+ 
A/i) where {AHq, A/i) are also chosen randomly. An exam- 
ple of computation is given in the next table. 



Iteration 


logiod^o -^ollw) 


logiodlM" "/^llw) 


1 


- .579029 


-1.358376 


2 


-3.003599 


-2.865026 


3 


-4.339497 


-4.122528 


4 


-8.234980 


-8.179398 


5 


-13.963299 


-14.029020 


6 


-14.022486 


-14.131066 



Here, we refind a pair {Ho,fi) starting from a 10% 
random perturbation. We see that the numerical convergence 
is obtained after 6 iterations. Note also that the quadratic 
convergence is observed. 

B. Continuation method 

In a second test, we use the continuation method pre- 
sented in IIV-CI to tackle a problem where the algorithms of 
Section IIV-BI do not apply. Given a target Utarget obtained 
with the field e and a pair {Ho,ii) that is chosen randomly, 
we look for the operators Hq and /x' that solve the control 



problem associated to the field cos(3t) and the target Utarget- 
The direct use of the Newton method of Section IIV-BI does 
not work: in this case, the algorithm does not converge. The 
continuation method enables to solve this problem. Using 
50 = 1/4, and 10 iterations of the Newton method as inner 
loop, a relevant pair {H^,^') is obtained. 
This example has been reproduced for numerous random 
initial pairs {Ho,fi). 
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